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Previously we identified the receding flow, where two fluid streams recede from each other, as an open 
numerical problem,’ because all well-known numerical fluxes give an anomalous temperature rise, thus called 
the overheating problem. This phenomenon, although presented in several textbooks,”* and many previous 
publications, has scarcely been satisfactorily addressed and the root cause of the overheating problem not well 
understood. We found that this temperature rise was solely connected to entropy rise and proposed to use 
the method of characteristics to eradicate the problem.”° However, the root cause of the entropy production 
was still unclear. In the present study, we identify the cause of this problem: the entropy rise is rooted in the 
pressure flux in a finite volume formulation and is implanted at the first time step. It is found theoretically 
inevitable for all existing numerical flux schemes used in the finite volume setting, as confirmed by numerical 
tests. This difficulty cannot be eliminated by manipulating time step, grid size, spatial accuracy, etc, although 
the rate of overheating depends on the the flux scheme used. Finally, we incorporate the entropy transport 
equation, in place of the energy equation, to ensure preservation of entropy, thus correcting this temperature 
anomaly. Its applicability is demonstrated for some relevant 1D and 2D problems. Thus, the present study 
validates that the entropy generated ab initio is the genesis of the overheating problem. 


I. Introduction 


HE overheating problem is one of the open CFD problems that yet to be understood and correctly computed; it has 

been found in two scenarios occurring in two opposite situations—one related to colliding compressive flows and the 
other to receding expansive flows. The former results either from reflection of a flow from a wall or equivalently from 
two colliding flows, manifested as a spike in temperature, and is first studied extensively by Noh.° This overheating 
problem still challenges CFD researchers today, see for example the books by Toro” and LeVeque.* The receding 
flow problem, while less known, occurs during rarefaction: for example, two streams recede from each other, thereby 
creating a rarefying region between them, where all thermodynamic properties—pressure, density, and temperature 
should be decreasing. Unexpectedly, a numerical anomaly arises for an isentropic expansion process during which the 
calculated temperature not only departs significantly from the correct profile, but a noticeably odd increase appears in 
the center region, as depicted in Fig. 1. Two typical solutions are displayed—between the diffusive HLLE method and 
the sharp Godunov method,” both showing the same oddity in temperature, although one is more smeared than the 
other. Toro” states ”... both pressure and density are close to zero and thus small errors will be exaggerated by their 
ratio” (p. 227). However, we showed in our previous studies'**> that this anomaly occurs even in early times when 
both pressure and density are still O(1). 

Based on the previous study,” the temperature rise is established to have a strong connection to the entropy rise. 
Hence, an alternative path is to study the entropy evolution for gaining insight into the problem. This motivates us 
to pursue the present study with the objective to analytically uncover the root cause why the entropy is generated by 
numerical fluxes in an isentropic flow. Then, an entropy-preserving system, which replaces the energy equation with 
the entropy equation, is implemented to confirm that the temperature anomaly is indeed eradicated. 

The paper is organized as follows. We first give a brief review of our previous attempts and findings, then we 
develop an analytical framework to look into the entropy change. Finally we incorporate the entropy equation in the 
solution procedure and confirm that if the entropy is preserved then the temperature anomaly disappears. 


*Senior Technologist, Aeropropulsion Division, MS 5-11, 21000 Brookpark Road, Cleveland, OH 44135. Fellow. 
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Figure 1. Receding flows computed on a grid of 100 cells, with the initial condition, (p, p, wu), = (1.0, 0.4, —2.0) and (p, p, vu) = (1.0, 0.4, 2.0) Ref. [°]. 


II. Statement of the Problem 


II.A. Definition of the overheating problem 


For the overheating problem discussed in this paper, we shall focus on the receding problem, specifically a Riemann 
problem with the following initial condition: 


UL < 0, PL, PL (1) 
ur > 0,pR,PR (2) 


where the “L” and “R” states may be of the same or different magnitudes. We shall confine our discussion to the case 
of equal magnitudes, uw, = —ur, pL = PR, pL = pr. The strength of rarefaction is defined by the relatively Mach 
number |Mp— M_|. 

Our interest in studying this seemly simple problem is motivated by the desire to understand the numerical mech- 
anism responsible for the overheating. A brief summary of our previous efforts in this regard will provide a useful 
guide for gaining insight into the problem. 


IL.B. Previous findings'*> 
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grid perturbed by 10%. 


Figure 2. A plane shock wave traveling in a constant-area channel in which the center grid line is perturbed at alternating even and odd points with a magnitude 
of 1% of grid spacing AY, results are obtained using the Godunov, Roe, and AUSM* -up with various CFL numbers. The AUSM*-up solution remains clean 
even with a 10% grid perturbation (Ref. Up. 


In the previous study of some ”difficult’ numerical problems, we characterize the difficulties into two categories: 
(1) fundamental and (2) operational. The ’operational” category refers to those problems remediable by change of 
operational procedures such as methods, mesh distribution, time step, limiters, etc. The carbuncle problem, which 
has remained a controversial topic, is an example of such operational problem, based on the evidence in our study: 
sample results shown in Fig.2 demonstrate that the shock anomalies associated with different methods can be made to 
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disappear by simply changing the time step. The difficulty in the ’fundamental” category, on the other hand, cannot 
be overcome by simply changing numerical parameters and/or switching numerical schemes. The receding problem 
considered in this paper belongs to this category. Figure 3 shows the typical calculated temperature profiles for 
subsonic and supersonic flows, in which again the temperature rises to a peak value at the center between two streams. 
From our previous studies, this phenomenon is universal irrespective of flux functions used, although quantitative 


differences exist among them. 
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(b) M=0.8, 10,000 cells. 
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(c) M=2.5, 100 cells. 


(d) M=2.5, 10,000 cells. 


Figure 3. Solutions of receding flows moving at M=0.8 and 2.5 by the Godunov method on 100 and 10,000 cells respectively (Ref. [']). 


Figure 4. Evolution of temperature and entropy at various time slices (Ref. [°]). 
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The incorrect temperature behavior is found to be related to entropy production,’ as depicted in Fig. 4. We see 
that (1) the entropy is increased at the first moment at the center, clearly violating the second law of thermodynamics 
for this smooth flow, and (2) this entropy violation continues to spread spatially. Furthermore, the entropy peak at 
the center reaches a maximum value and then stays at a saturated value, which is found to be rather insensitive to the 
number of mesh points used. This explains why the overheating is not shrunk by changing the order of accuracy or 
mesh size. It also strongly suggests that in order to resolve this unphysical overheating problem in the receding flow, 
preservation of the entropy is the key. This viewpoint is further confirmed that if the entropy condition is enforced, 
then the temperature indeed behaves correctly, as will be shown later. 

The link of the temperature increase to the entropy production is further analyzed subsequently,’ in which it is 
argued that the overheating is inevitable in the current finite volume discretization framework. This link of entropy 
production and temperature rise is found to follow a numerical heat diffusion model, 


Tds = kdT (3) 


T(a+ Az) 
T(x) 
Figure 5 shows the variation of these two calculated quantities between two neighboring points for the entire 

computation domain, confirming that the relationship in Eq. 4 holds, the entropy increases with temperature, and the 

slope manifested as the numerical diffusivity of this calculation is nearly a constant. The points further away from the 
center correspond to the region of large changes in entropy and temperature, i.e., those near the center of the receding 
flow. 


As = s(a + Az) — s(x) = Klog (4) 
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Figure 5. Temperature variation A,, log T vs. entropy variation A, s. 


To remove this temperature anomaly in the receding flow, the method of characteristics (MOC) was used to replace 
the conventional discretized system to advance the solution; this approach has proven to be successful in correcting 


the temperature overheating.” 
However, the reason why the entropy is produced in the first place remains a mystery. This motivates the current 


study, as will be described in what follows. 


Ill. Formulation 


As evident from the previous study, the entropy is generated and plays a unique role in the overheating problem. 
The key is to find how the entropy is produced in the discretized conservation laws employed for the numerical solution. 
A fundamental thermodynamic equation for an ideal gas gives 


d d Ve d 
Tds = de — 2dp =e(2 - 7) = OF ap Cur dies) (5) 
p Pp p Pp p p? 
Or C 
Pp v 2 
ds = Cy d(log —) = —* (dp — a?d (6) 
( 8 oy? Pa p) 
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The first equality in Eq. 6 allows us to define the usual entropy variable for an ideal gas as, (omitting an integration 
constant without loss of generality) 
R 
s=Oylog 2, C, = —— (7) 
p? hee 
where R is the gas constant. The second equality in Eq. 6 allows to relate the time rate of change of entropy to that of 
the conservation laws via pressure and density. From the energy equation, we find 


OpE  O (pu)*) it ie? 290) -( i MEO Opu (8) 
ae Ot Op 7 1 ot” Ot’ Sy 12 a 
Combining Eqs. 6 and 8 yields 
pds OpE a? u*. Op Opu 
= ( ) u (9) 
Rot ot ‘y-1 2’dat at 


Clearly this equation characterizes the balance of entropy change and the changes of mass, momentum, and total 
energy. This is a very illuminating equation, from which we see how the entropy can be produced. Let us, for 
the purpose of illustration, consider the 1D conservation laws. The conclusion reached herein however holds for a 
multidimensional system, as will be shown later with some numerical results. 


Op Opu 
a On 
Opus 0, g 
OE Bg ee +p) (11) 
OpE OpuHl 
pote ae 12 
Ot Ox - 


Now substituting the above spatial derivatives in the balance equation Eq. 9, the continuous system becomes, after a 
simple algebraic manipulation, 


2S gy 1 
Rat Rk dx a 
Since p/R # 0, it simplifies to 
Os Os 
FE + us = 0 (14) 


This is the familiar entropy equation, leading to the well-known result for a smooth flow—The entropy remains un- 
changed along the particle path. 

The receding flow considered above is an isentropic flow initially and it should theoretically remain isentropic 
throughout according to Eq. 14. However all the numerical solutions predicted otherwise. The key difference between 
the theoretical and numerical results lies in the fact that the discrete version of Eqs. 10-12 is not consistent with the 
continuous system, as will become evident below. 

Let us consider the finite volume discretization of Eq. 9, in which the right hand side is replaced by numerical 
fluxes at the cell faces, see the following schematic. 


ion 


-1/2 1/2 3/2 


Figure 6. Schematic of finite-volume cells and indices. 
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For the “R” cell bounded by faces “3/2” and “1/2”, with initial conditions: up = —uz, = uo > 0,prR = pt = 
fo, PR = PL = Po, we obtain 
p Os a2 ua 


= Roe = —[(eu)sj2— (oul )i2)+(— — 17 Meu)sy2—(eu)aa] +ug[(pu?+p)3/2—(pu+p)1/2] (15) 


where the coefficients associated with the momentum and mass flux differences are evaluated at the initial time. The 
fluxes for the face “3/2” is simple, based on the consistency requirement (F' (Ur, Ur) = F(Up)), we have 


a U 
(puH)3/2 = pototo = pouol— . it a 
(16) 
(pu)3/2 = pouvo 
(pu? + p)s/2 = pod + Po 
And for the face “1/2” where we have u 1/2 = 0, by virtue of anti-symmetry, the fluxes are 

(puf1)1/2 = 0 

(pu)i/2 = 0 (17) 


(pu? + p)ij2 = ™M4/2 + Pie 


where the momentum flux m,/ and pressure flux p,/2 can take different forms for different flux schemes; they are 
given in Table 1. For all the numerical fluxes tested, we find 


Pij2=€po, 1>e=0 (18) 


Note that the € for the exact Riemann solver is valid up to the vacuum condition (Mp = 2/(7 — 1) in this case). The 
values of € of some flux schemes for Mp = 0.8 are given in Table | for reference. By substitution of Eqs. 16 and 
eq:1/2, Eq. 15 reduces to 


Roe = uo|(po — P1/2) — M1/2] = uol[(1 — €)po — ™y/2] (19) 
“Rn 

The first term is non-negative, thus contributing an increase in entropy (since p/R > 0), while the second term may 
lead to a decrease in entropy if m2 > 0. This latter situation happens only when ug > Go for the Roe scheme, but this 
results in violation of the entropy condition, which in turn leads to a negative temperature and failure of computation. 
Hence, this entropy (thus temperature) positivity condition requires that ™m 1/2 < 0. Notice that the van Leer flux gives 
my1/2 < 0, while all other schemes, as shown in Table 1, give m1/2 = 0,VMo. Hence, for an entropy-satisfying 
(As > 0) flux scheme (implying m1/2 < 0), the entropy is generated within the first time step in cell ’R”. Similarly 
we get the identical result for the ’L” cell bounded by faces ”1/2” and ”-1/2”, 


O 
¢ Bae = Moll ~ Po ~ maya (20) 


with 


(puH)_1/2 = pouo( 


(pu)—1/2 = —pouto eh 


(pu? + p)-1/2 pou + Po 


This shows that the entropy increase occurs in both the ”L” and ”R” cells, symmetrically with respect to the face 
”1/2” straddling them. Based on the analysis developed above, it is the numerical flux, specifically the pressure flux, 
at the cell face that creates the entropy. The numerical pressure flux at the interface ”1/2” results in an imbalance of 
pressures on the opposite faces, thus creating a net increase in entropy. This imbalance in pressure will continue to 
develop in the neighboring cells at the following time steps so long as the flow is evolving (changing) spatially. Hence, 
the entropy generation and temperature rise continue to spread spatialwise. 
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Table 1. Expression of m1 /2 and p; /2 (or €) of « in Eq. 18 for various numerical fluxes, which can be found in Ref. Pe 


Flux Scheme Exact Riemann AUSM-family Roe Van Leer | Steger-Warming 
M41 /2 0 ; 0 potio(uo — @o) | Mi2,v4b 0 
€,Pij2 = Po | (L— 23*Mo)7-! | €1/2,ausm+ 1 €1/2,VL €1/2,S—W 
(Mo = 0.8) 0.295 0.071 1 0.056 0.040 
1 2 sey, 
—5pougao(Mo = 1) 4 if My < L, 
m = 4 22 
ewe 0, otherwise. a 
2(Mop —1)?(2+ Mo), if Mo <1, 
ages at eo) NO (23) 
0, otherwise. 
1(Mp —1)7(2+M, 3Mo(M? —1)?, if Mo <1, 
ere are 3 (Mo — 1)?( 0) — §Mo(Mo — 1) 0< (24) 
0, otherwise. 
(My —1)?, if Mp <1, 
€ “w= 25 
ener 0, otherwise. oy 


From the above result, we now conclude the following. 


Theorem: Under a finite volume discretization, the entropy production for a receding flow has its origin in the 
pressure flux, established ab initio within the first time step. This is true regardless of cell size, flux function scheme, 
or time step; however the amount of entropy produced will depend on the flux functions used and the receding velocity. 

Remark 1: For an entropy-satisfying condition, i.e., 0s/Ot > 0, the momentum flux m 1,2 must be non-positive; 
otherwise a negative temperature will result, rendering the calculation terminated. This happens to the Roe fluxe in 
supersonic receding flow. 

Remark 2: The above analysis can serve as a simple tool to examine whether a flux scheme satisfies the entropy 
or positivity condition. 


Although the above analysis only proves that the entropy is generated within the first time step, our calculations 
confirm that the trend continues in time and extends spatially. In fact, we find from the calculations in ['] that for 
t > 0: (1) the overheating error remains even after increasing the number of cells by 100 times and (2) the error 
increases with the magnitude of the separation velocity (the strength of rarefaction), also evident in Fig. 3. 

We re-iterate that the continuous system does guarantee preservation of entropy, it is only in the spatially discretized 
system that the entropy is produced. With the root cause established, the remaining question as to how to eradicate 
this numerical anomaly becomes clear. 

Clearly, in order to eliminate this entropy production and the overheating anomaly the usual finite volume frame- 
work has to be revisited and modified. An alternative is to use an approach that obeys the entropy equation. The 
method of characteristics (MOC) is one that inherently includes the entropy condition. The MOC was shown to be 
successful in eradicating the overheating problem in some relevant flow problems.° However, the MOC is rather cum- 
bersome in dealing with geometry and boundary conditions; moreover it is hard to extend to higher order (temporal 
and spatial ) accuracy. 

If all is possible, we still prefer the conventional finite-volume formulation, the only requirement is to incorporate 
the entropy condition. This can be achieved by using the entropy equation to replace the energy equation. 


TiI.A. Use of the entropy equation 


The nonconservative entropy equation*, Eq. 14, can be used along with the continuity and momentum equations in 
conservation form, and a simple upwind discretization is used. 


_ At 
Ar 


“Since an isentropic flow implies a smooth flow, conservation form is not necessary. 


n+1_ on 
i 33 
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(26) 
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Or a conservative form can be used by combining with the continuity equation, given as:° 


Ops Opus 

at Om 
It now takes the place of the energy equation; it is simple and cast in the same framework as the other conservation 
equations. For example, the entropy flux (pus) at the cell interface is evaluated exactly the same way as the other 
fluxes for the AUSM-family schemes. Figure 7 demonstrates that the result of enforcing the entropy equation for the 
1D receding flow—maintaining constant entropy and eliminating the temperature rise; similar results are also obtained 
for other flux schemes. We also note that the conservative form of the entropy equation gives results indistinguishable 
from that of the nonconservative form. 


=0 (27) 


R-/-R R-/-R 
100 cells 100 cells 

AUSM*-up CFL=0.50 N= 60 AUSM*-up CFL=0.50 N= 60 

1.5 —4— Pressure 1.5} —4— Pressure 

—— Density ——Density 
—=—Temperature —s—Temperature 
—— Mach (rescaled) —— Mach (rescaled) 
—— Entropy —— Entropy 


(a) AUSM*+-up (b) AUSM*+-up with s-equation 


Figure 7. Effects of enforcing the entropy equation. (a) AUSM* -up and (b) AUSM* -up with the entropy equation for the 1D receding flow problem (M=0.8). 
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Figure 8. Effects of enforcing the entropy equation. (a) AUSM* -up and (b) AUSM* -up with the entropy equation for a 2D receding flow problem (M=0.8). 


>The conservation form is not unique as it can in general appear as Ops™ /Ot + Opus™ /Ox = 0, m = 1,2,... 
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In what follows we shall consider some relevant 2D rarefaction problems. The 2D version of the entropy equation 
is applied. First we calculate a 2D receding flow, defined by the initial condition: 


UR up #0, vR=—v_,=0, pR=PL, PR=PL; (28) 


where wu and v are the x and y components of the velocity, respectively. The flow is calculated on a 2D slanted 
mesh. The results obtained by using the energy equation and the entropy equation are given in Fig. 8. As in previous 
examples, overheating is seen with the conventional formulation, while the entropy equation gives rise to the expected 
correct result in which the entropy remains constant and the temperature behaves properly. 
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Figure 9. 2D Prandtl-Meyer flow (M=2.0, corner angle of 20°. (a) AUSM*-up (100 x 400 mesh) and (b) AUSM*-up with the entropy equation (100 x 100 
mesh); reference solution is obtained with a 100 x 400 mesh. 


Next we consider a supersonic (M=2.0) flow over an expansion corner of 20°, or the Prandtl-Meyer flow. Here 
we show the calculated results without and with using the entropy equation; the inaccurate prediction of the wall 
temperature, namely overheating, is clearly seen in Fig. 9 when the conventional conservation equations are employed, 
even with a very fine mesh of 100 x 400. On the other hand, the entropy equation gives a correct wall temperature and 
excellent agreement with the reference solution, with a coarse mesh of 100 x 100. The reference solution refers to the 
isentropic solution on the mesh of 100 x 400. 

Table 2 gives the comparison of calculated wall temperature with different grid sizes; it shows that the result is 
hardly improved by reducing the grid spacing by eight times when the entropy condition is violated; the error from the 
reference solution remains about 10%. Interestingly, the isentropic solution on the coarsest mesh (100) already gives 
a wall temperature essentially identical to the reference value. 


Table 2. Comparison of predicted wall temperature with different grid sizes to the reference value. 


Grid points | N, =50 | N, = 100 | N, = 200 | N, = 400 | s-eqn, N, = 50 | Reference solution 
Twall 0.7793 0.7767 0.7754 0.7742 0.6893 0.6900 


It must be reminded that while the use of the entropy equation cures the overheating problems, its applicability 
is limited to smooth flows only. Even within the realm of inviscid flow, many are nonsmooth and involve shocks. 
Including the entropy equation in the system will still yield discontinuous solutions across a shock but with entropy 
unchanged, hence its strength is incorrect. For transonic flows, the error is small, (O(/ — 1)3), this formulation 
nevertheless is not general. We use the entropy equation only to provide the numerical evidence for the analysis 
developed above that the entropy generation is the reason for the temperature anomaly. 
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IV. Conclusion 


We have achieved the objective of this study: to uncover the root cause of the overheating in the smooth receding 
problem. The root cause of the overheating has its origin in the entropy generation ab initio, which is tied to the 
pressure flux in the finite volume discretization of the conservation laws. The overheating is found to be inevitable for 
all numerical flux schemes used, as confirmed in the numerical tests. The rate of overheating is proportional to the 
flow separation speed and is a function of flux scheme. 

To support the conclusion of the analysis, an isentropic system is used: the entropy equation is employed to replace 
the energy equation. The use of the entropy equation proves that incorporating the entropy condition in the discretized 
system eradicates the temperature overheating problem, ensuring the correct temperature and entropy distributions. Its 
efficacy has been demonstrated in several relevant rarefaction problems. Especially the inaccuracy in the calculated 
wall temperature for the 2D Prandtl-Meyer problem by the conventional finite volume approach may have a practical 
implication. 

Finally, we also remind that the applicability of the entropy equation is limited to isentropic flows which in practice 
are rather rare. Its usage in this study is only to provide numerical evidence of the root cause of the overheating problem 
and is not recommended for general problems. However, if an occasion of inaccurate temperature prediction arises for 
an isentropic flow, we now know the cause and how to fix it. 
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